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We consider (2+l)-dimensional beams, whose transverse size may be comparable to or smaller 
than the carrier wavelength, on the basis of an extended version of the nonlinear Schrodinger equa- 
tion derived from the Maxwell's equations. As this equation is very cumbersome, we also study, in 
parallel to it, its simplified version which keeps the most essential term: the term which accounts 
for the nonlinear diffraction. The full equation additionally includes terms generated by a devia- 
tion from the paraxial approximation and by a longitudinal electric-field component in the beam. 
' Solitary-wave stationary solutions to both the full and simplified equations are found, treating the 

Qh i terms which modify the nonlinear Schrodinger equation as perturbations. Within the framework 

pH . of the perturbative approach, a conserved power of the beam is obtained in an explicit form. It is 

found that the nonlinear diffraction affects stationary beams much stronger than nonparaxiality and 
longitudinal field. Stability of the beams is directly tested by simulating the simplified equation, 
with initial configurations taken as predicted by the perturbation theory. The numerically generated 
solitary beams are always stable and never start to collapse, although they display periodic internal 
vibrations, whose amplitude decreases with the increase of the beam power. 
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I. INTRODUCTION 



Self-focusing of optical beams in nonlinear media |^ is a subject of intensive study since the critical nature of this 
^ , phenomenon (a collapse of the beam) has been identified in the first studies (|-|^. Different theoretical approaches 
have been developed for the analysis of self- focusing, including the moment theory J^, variational method (see, e.g., 
and references therein), paraxial-ray approximation "adiabatic" description |lO), etc. With a probability for a 
collapse, no stable self-trapped beams can exist in bulk Kerr-mcdia. Therefore, a considerable part of the studies in 
5— i the field has been directed towards finding mechanisms that limit the development of the collapse and thus provide 
, for the existence of stable self-trapped beams in two transverse dimensions. It has been shown that the saturation 
of the nonlinearity of the medium arrests the collapse and stable two-dimensional self-trapped beams can indeed exist 
in such media The plasma generation by a self- focusing laser beam fl^-p^ through multiphoton and avalanche 
ionization, with its negative (exponential) contribution to the refractive index of the media, can also stop the collapse 
by leading to formation of filaments. 

However, as it is pointed out in Ref. the occurrence of the self-focusing in different media (gases, liquids, 
solids) makes it insufficient to identify a stabilizing effect for each separate case, and calls to search for universal, i.e., 
medium-independent, mechanisms which yield nonsingular behavior of the beams. In this connection, it is necessary 
to stress on that the dynamics of narrow beams described by the full system of the Maxwell's equations may be 
drastically different from the predictions of the usual parabolic nonlinear Schrodinger (NLS) approximation valid for 
broad beams (see, e.g., jl^ ] and references therein). In fact, the Maxwell's equations implicitly contain few different 
mechanisms that can check the collapse of the narrow beams. 

One of such mechanisms is nonparaxiality of the beam pro pag ation, which becomes important at the advanced 
stages of the self- focusing. It was studied both numerically Pq,[17|-[l9|] and analytically pO], concluding that the 
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nonparaxiality replaces the catastrophic focusing with a sequence of focusing-defocusing cycles. Beam filamentation 
due to the nonparaxiality was demonstrated as well jl^ . 

It was also shown that when the vector nature of the beam field is taken into regard along with the nonparaxiality 
the peak intensity occuring in the self- focusing is almost ten times lower than that produced by the nonparax- 
iality alone for the same input beam, suggesting that other effects, which limit the self-focussing stronger than the 
nonparaxiality, could exist. Recent analysis ||2^,^ has demonstrated that the efficiency of the vector model in limiting 
the narrowing of the beam in the (l-l-l)-dimensional case is due not to the small longitudinal electric- field component 
but rather to a scalar term which is present in the model. This scalar term, accounting for the rate of the spatial 
variation of the nonlinear polarization, combines effects of nonlinearity and diffraction. A notation "nonlinearly in- 
duced diffraction" was introduced 23| to stress on this factor which restricts the beam narrowing predicted by 
the usual NLS equation. In this connection it should be mentioned that a fundamental limit for the soliton width 
has been found [2^p6| to exist in the " subwavelength" range (beams with a transverse size comparable with or 
smaller than the carrier wavelength) due to the terms that come from Maxwell's equations but are absent in the usual 



(l-l-l)-dimensional NLS equation. Given the results in Ref. |2ll and conclusions drawn in Refs. |22 2q], the existence 



of a very narrow stable self-trapped beam in an ideal bulk Kerr-medium is still an open problem which needs to be 
addressed, which is the subject of the present work. 

The rest of the paper is organized as follows. In Section 2, we derive, from the Maxwell's equations, an equation 
to govern the evolution of narrow (2-|-l)-dimensional beams in Kerr media. This is a very cumbersome generalization 
of the NLS equation, with extra terms representing three different effects: nonlinear diffraction, nonparaxiality, and 
the influence of the small longitudinal component of the electric field. In view of the very complex character of the 
equation, we also consider a simplified version, that keeps only the nonlinear diffraction as a new effect. In section 3, 
we develop a perturbation theory which treats all the new terms in both equations as a perturbation, which is true in 
the case when the beam is still relatively broad. Further analysis in section 4 shows that the effect produced by the 
nonlinear diffraction is much stronger than those generated by the nonparaxiality and longitudinal field, which justifies 
the use of the simplified equation, at least for beams that are not extremely narrow. Another noteworthy result of 
the perturbation theory is that, although the terms which modify the NLS equation destroy the circular symmetry 
of the beam, the deviation from the symmetry is, in fact, fairly weak. In section 5, we show that the perturbative 
approach also allows one to find an expression for the conserved power of the beam, which is not a trivial issue in 
the present model. Lastly, in section 6 we display result of direct numerical simulations of the beam propagation, 
starting from initial configurations predicted by the perturbation theory. The direct simulations are performed within 
the framework of the simplified equation only. As a result, we conclude that the beams never show collapse and are 
always stable, showing small internal vibrations, which are generated by a deviation of the initial configuration from 
an exact solitary-beam shape. 



II. THE BEAM PROPAGATION EQUATION 



We consider the propagation - along the z-axis - of a light-beam, extended in the {x, y)-plane, in a bulk isotropic 
Kerr-medium, characterized by the usual relation between the nonlinear polarization P^^ and electric filed E, P^^ — 
Eqx ij^x + ^z) where x is the Xxxxx — Xzzzz component of the third-order susceptibility tensor [ p7| . The wave 
equation, 
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yield a scalar nonlinear wave equation in the following form: 
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In Eq.(3), lo is the frequency of the carrier wave, E = E^, A± = + {d"^ /dy"^) is the transversal Laplace 

operator, el and [3 = {lo /c) £l are the linear permittivity and the linear propagation constant (wavenumber) , respec- 
tively. After rescahng 2[3x x, 2f3y y, 2(3z z, jE E, with 7^ = (Sw^x) / (l6/?^c^), and some manipulation, 
Eq. (3) transforms into 
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The modifications of the (2-|-l)-dimensional NLS equation involved in Eq. (Q) stem from nonparaxiality (the fourth 
term) and nonlinear diffraction and longitudinal field component included in a combined manner in the last three 
terms. It makes sense to consider Eq. (0) parallel to its "model," simplified, version 
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where only the nonlinear diffraction [the last term in Eq. (H)] is taken into account to modify the NLS equation 
The relevance of approximating Eq. (M) by Eq. (H) will be verified a posteriori. 



III. SOLITARY- WAVE SOLUTIONS AND PERTURBATION THEORY 



Seeking for solitary-wave solutions to Eqs. (Q) and (|^) in the form E{x,y,z) ~ F{x,y)exp{ifiz), where /i is a 
nonlinear wavenumber shift of the solitary-wave solutions considered, one obtains 
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In the case of Eq. (|^), a = 12, 6 = 24, and k = fi, while Eq. (|^) corresponds to a = 32/3, b = 68/3 and k = fi + . 
Besides the linear diffraction in the second {y-) transverse direction, Eq. exactly coincides with the one-dimensional 
equation derived in Ref. p^for both the vector nonparaxial model and for the scalar "nonlinear diffraction" model. 

Generally speaking, Eq. (^) must be solved numerically as it stands. However, for the case of small gradients a 
perturbation scheme can be developed, which simplifies the problem considerably. In fact, this implies seeking for 
first corrections to the usual broad solitons, generated by the new terms obtained from the Maxwell's equations. To 
this end, F{x, y) and k are presented in the form F{x^ y) = Fq{x, y) -\- f{x, y), k = kg + Ak, where / and Afc are the 
first order corrections, so that the substitution of this into Eq. (^) yields 
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In Eq. (|^), ko — jiQ is the zeroth-order approximation for the nonlinear wavenumber shift. The transformation to the 
polar coordinates x = rcos(p, y — rsimp transforms Eqs. and (^) into 
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respectively. Equation (^), which is azimuthally symmetric, gives rise to the classical Townes soliton while Eq. 
(10) generates a first-order correction to it. In Eq. (10) the variables can be separated by means of the substitution 
/(r, (fi) — fo{r) + fi{r) cos 2ip, the resulting equations for /o and /i being 

dr2 + r dr 1" ~ 2 [ dr^ r dr 2 [ dr 



Thus, the first-order perturbative solution to Eq. (pf) amounts to solving three ordinary differential equations, viz.. 



Eqs. (|9|), ( |ll| ) and ([l2|). This can be easily done by means of the shooting method |28|. We supplement the equations 
with obvious boundary conditions that single out soliton solutions: 

dFojr^O) ^ dfojr = 0) ^ dhjr ^ 0) _^ 
dr dr dr ' 

Foir ^0)=A> 0,/o(r - 0) - /i(r - 0) - 0, 

Fo{r ^ w) ^ 0, /o (r ^ oo) ^ 0, /i (r ^ oo) ^ 0. 



IV. RESULTS OF THE PERTURBATION THEORY 



The solution of Eq. (9) is the same as that obtained in Ref. For this solution, fep ~ and the power 

oo 

P that it carries is P = 27r J rF^dr « 11.70, which is the known "critical power" for the weak collapse in the 

two-dimensional NLS equation [PIItH^. 

Solutions to Eqs. ( pi] ) and ( p^ are displayed in Figs. |l](a) and|l](b), respectively. The solutions marked by "NLD" 
pertain to the case of a = 12, 6 = 24, which corresponds to the simplified equation (^, and those marked by "vector" 
pertain to a = 32/3 ~ 10.67, h = 68/3 ~ 22.67, which correspond to the full equation (^. For comparison, the 
zcroth-order solution (the Townes soliton) is shown in Fig. [|(c). It is obvious that the nonlinear diffraction is the 
main factor which causes deviations from the solution of the NLS equation and the effect of the longitudinal field 
component is quite small. 

In Fig. |2[ results for the half- widths Rx , Ry of the beam (defined at half-maximum of its amplitude) are shown. 
It can again be concluded that the nonlinear diffraction term is the main reason which restricts the narrowing of 
the Townes soliton with the increase of its on-axis intensity (denoted as "standard" in Fig. ||). The influence of 
the longitudinal field component shows up into a small reduction of the effect of the nonlinear diffraction. It is also 
worth noting that the half-width Rx in the direction of the main electric-filed component (x) is much more affected 
by the nonlinear-diffraction term than that {Ry) in the other transverse direction (y). The latter feature characterizes 
an effective anisotropy induced by the new terms. Note that deviations from the circular symmetry of a beam 
propagating according to an equation similar to Eq. (^ was reported in Refs. [p]pl|]. 

The comparison of the contributions to the nonlinear wavenumber shift fi stemming from nonparaxiality, nonlinear 
diffraction and longitudinal electric field offers an additional way to estimate a relative strength of effects generated 
by these factors. In the case when the simplified model (^ is used, which keeps the nonlinear diffraction but 
neglects the nonparaxiality and the influence of the longitudinal field, i.e., a = 12 and 6 = 24 in Eq. (§), the 
result is fiNLiD ~ + Yn *° A^-order. On the other hand, the full model corresponding to Eq. (||), 

i.e., with a — 32/3 and b = 68/3 in Eq. (^, yields fivN ~ + Y74 (l — 20^52 ~ 13 57 ) • ^^^^ correction in 

the parenthesis is due to the vectorial structure of the field, and the second one is an effect of the nonparaxiality. 
This result demonstrates, that within the framework of the first-order perturbation theory, the contribution to the 
nonlinear wavenumber shift stemming from the nonlinear diffraction exceeds those originating from nonparaxiality 
and longitudinal field component by approximately 13 and 20 times, respectively, thus validating the use of Eq. (|^) 
instead of the full equation (^) for the study of the solitary waves. 
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V. AN APPROXIMATE POWER-CONSERVATION LAW 



According to the results of the perturbation theory for the case of small gradients (Figs. and g), the deviation 
of the shape of the beam from the circular symmetry is small (although existing) . Introducing the polar coordinates 
into Eq. (^ and assuming weak dependence of the field envelop on the angular variable Lp, one obtains 

Multiplying Eq. ( p^ ) by (^E* + \E\'^E*^, subtracting the complex conjugate and integrating over the transversal plane 
while keeping the first order corrections only, a power conservation law in the form 

27r 00 


is obtained. Note that the second term in the integrand is due to the nonlinear diffraction which comes to involve the 
power density associated with the nonlinear polarization into the power conservation law. 

In Fig. ^ the dependence of the power invariant, given by Eq. (14), on the on-axis intensity of the beam is shown. 

For comparison, the same dependence for the case of the standard NLS equation (/ / r\E\'^drdip) is also presented. 



It is evident that the solitary waves studied here exist for power levels above the critical. The fact that the power is 
a growing function of the on-axis intensity suggests stability for the solitary wave solutions (see also the discussions 
m Ref. |l| in the case of the scalar nonparaxial model). 



VI. NUMERICAL SIMULATIONS OF THE BEAM PROPAGATION 



As Fig. 1^ shows, the presence of nonlinear diffraction increases the power necessary for self-trapping (i.e., the power 
at which the diffraction and nonlinearity balance each other). This is demonstrated by the numerical simulations 
displayed in Fig. |^(a), where an initial configuration in the form of the Townes soliton, see Eq. (|9|), propagates 
according to the Eq. (g). We do not simulate the full equation (^), which has an extremely complex form, since the 
results presented above strongly suggest that the simplified equation (^ may describe essential features of the beam 
dynamics quite adequately. Figure H shows that in this case the beam diffracts away, while, as it is well known, the 
usual NLS equation predicts collapse for the same initial conditions [see Fig. ^(c)]. In Fig. ^(b), the approximate 

00 oc 

invariant given by Eq. (|lj) and the invariant of the NLS equation in its standard form, J J \E\'^dxdy^ are shown 

— 00 — oc 

vrs. the propagation distance for the case when Eq. (||) is solved. The fact that the NLS invariant increases when 
the beam diffracts can be explained in the following way. As it has been mentioned above, the full invariant, given 
by Eq. (^4|), takes into account the power density associated with both the linear and nonlinear polarization of the 
medium and is conserved^ provided the gradients remain small enough. When the beam diffracts, the contribution of 
the nonlinear polarization decreases, transferring the power it carries to the linear-polarization part, which therefore 
increases in accord with what is seen in Fig. ^(b). When the beam is self- focusing, the computations show that the 
NLS invariant decreases, as one should expect. 

In order to test directly the dynamical properties of the solitary waves in the present model (first of all, their stability 
against the collapse), equation (^ is solved numerically with initial conditions produced by the perturbation scheme 
developed above [i.e., the solitary-wave solutions generated by Eqs. (^ - ( p^ ) have been used as initial configurations]. 
Figure ^ shows results for the evolution of the beam intensity and of the approximate power invariant given by Eq. 
(|lj). For comparison, the usual NLS invariant is also displayed. As is seen in Fig. ^, the beams experience periodic 
focusing and defocusing with a relatively small amplitude of the corresponding internal vibrations, rather than the 
collapse that would occur in the case of the NLS equation. The period of the oscillations rapidly decreases with the 
increase of the input power. 

The results presented in Figs. 5(a) and 5(b) are in a qualitative agreement with those in Ref. [|l^ where the 
non-catastrophycal propagation of Gaussian beams within the scalar nonparaxial model has been obtained. The 
soliton-like propagation studied later [^ within the same model shows the same type of oscillating behaviour. The 
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influence of the deviations of the initial conditions from the exact sohton-hke solution has been also checked. It seems 
that the vibrations obtained within the nonparaxial model in Refs. ||l8| , |l9|] and here, in a model which stresses on the 
nonlinearly induced diffraction, are general behaviour which is a leftover of the former collapse instability. 

The oscillations observed here can be related to the fact that the initial conditions, predicted by the perturbation 
theory, generate a nearly-solitonic (but fairly stable) beam, rather than an exact stationary soliton (cf. a similar 
situation in the multidimensional model with the quadratic nonlinearity p9|] ). However as the power of the beam 
decreases and, accordingly, the beam radius increases, the accuracy of the perturbation theory improves, hence the 
discrepancy between the exact stationary self-trapped beams and those predicted by the perturbative approach should 
decrease, i.e., lower-power beams should oscillate with a smaller amplitude. This, however, is not the case. A closer 
look at Fig. || shows that when the power is increased [from (a) to (b)], the relative amplitude of the oscillations 
decreases slightly. To check this trend, some results for the beam propagation at lower power levels are presented in 
Fig. ^. The three different curves in this figure correspond to initial conditions in which, respectively, 6, 14, and 20 
digits after the decimal point are kept to accurately implement the prediction provided by the perturbation theory. 
The corresponding exact values of the power are also given for each case. It is well seen from Fig. 6 that the evolution 
of the beam is quite sensitive to fine details of the initial conditions. In order to achieve stationary propagation, the 
initial conditions must be set with a very high accuracy, especially for lower power beams. It turns out that the latter 
(wide) beams are more sensitive to small deviations of their initial shape from exact solitary waves than the higher 
power (narrow) beams. So, we infer that, quite naturally, the stabilizing effect provided by the nonlinear diffraction 
gets enhanced with the decrease of the beam size and the growth of its power. 



VII. CONCLUSIONS 



Propagation of (2+l)-dimensional narrow beams in Kerr self-focusing media is studied on the basis of an extended 
version of the NLS equation derived directly from the Maxwell's equations. As this equation is very cumbersome, we 
also studied, in parallel to it, its simplified version which keeps the most essential new term, viz., the one accounting 
for nonlinear diffraction. The full equation contains, besides the nonlinear diffraction term, terms generated by a 
deviation from the paraxial approximation and by a longitudinal electric-field component in the beam. Solitary-wave 
solutions to both the full and simplified equations are found, first, by means of a perturbation theory which assumes 
that the new terms are still small, provided that the beam is not too narrow. Within the framework of the perturbative 
approach, small (but explicit) deviations of the beam from circular symmetry are found, and the power-conservation 
law is derived in an explicit form. It is shown that the nonlinear diffraction affects the stationary beams much stronger 
than the nonparaxiality and longitudinal field. Stability of the beams is tested by direct simulations of the simplified 
equation, with initial configurations taken as predicted by the perturbation theory. Numerically generated solitary 
beams are found to be always stable and never sliding into collapse, although they display periodic internal vibrations, 
whose amplitude decreases with the increase of the beam power. 
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FIG. 1. Solutions of Eqs. ( p4| ) and (^[), respectively, in (a) and (b) with a = 32/3, b — 68/3 (notation "vector") and a = 12, 
6 = 24 (notation "NLD") as these quantities specify the difference between Eqs. (^) and (g). In (c), the solution of Eq. (^ is 
presented together with the solutions fo and /i from (a) and (b) taken for the case of a = 12, b — 24. 

FIG. 2. Dependence of the half widths Rx and Ry on the on-axis intensity [A = Fo{r = 0)] of the beam. The values of 
the parameters are the same as denoted in Fig. 1(a). 

FIG. 3. Dependence of the power invariant P on the on-axis intensity A as given by Eq. (Q [NLD model, Eq. (|)] and by 
the NLS equation. 

FIG. 4. On-axis intensity vrs. the propagation distance obtained as a solution of Eq. (^ in (a) and its corresponding power 
invariant in (b); the power invariant corresponding to the NLS equation is given in (b) in order to show its changes along the 
propagation distance. In (c), the corresponding solutiuons of the NLS equation are shown. A = 0.11 in the initial condition. 

FIG. 5. Evolution of the normalized on-axis intensity [in (a) and (b)] as solutions of Eq. (^) and of the power invariant [in 
(c) and (d)] according to Eq. (p^). For comparison the variations along the propagation distance of the power invariant of the 
NLS equation are shown in (c) and (d). 



FIG. 6. On-axis intensity vrs the propagation distance: check of sensitivity with respect to the initial conditions. 
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